R Notebook

library(Seurat)
library(Signac)
library(EnsDb.Mmusculus.v79)
library(ggplot2)

counts <- Read10X_h5(paste0(params$out_prefix,"scATAC_P50/data/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5"))
metadata <- read.csv(
  file = paste0(params$out_prefix,"scATAC_P50/data/atac_v1_adult_brain_fresh_5k_singlecell.csv"),
  header = TRUE,
  row.names = 1
)

fragment.path <- paste0(params$out_prefix,"scATAC_P50/data/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz")

brain <- CreateSeuratObject(
  counts = counts,
  assay = 'peaks',
  project = 'ATAC',
  min.cells = 1,
  meta.data = metadata
)

Some cells in meta.data not present in provided counts matrix.

brain <- SetFragments(
  object = brain,
  file = fragment.path
)

brain <- RunTFIDF(brain)

Performing TF-IDF normalization

brain <- FindTopFeatures(brain, min.cutoff = 'q0')
brain <- RunSVD(
  object = brain,
  assay = 'peaks',
  reduction.key = 'LSI_',
  reduction.name = 'lsi'
)

Running SVD
Scaling cell embeddings

brain <- RunUMAP(object = brain, reduction = 'lsi', dims = 2:30)

11:57:47 UMAP embedding parameters a = 0.9922 b = 1.112
11:57:47 Read 3880 rows and found 29 numeric columns
11:57:47 Using Annoy for neighbor search, n_neighbors = 30
11:57:47 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:57:47 Writing NN index file to temp file /var/folders/zc/29jxbf153g37b_cjp6_gkxnxxh967n/T//RtmpZ9ef1K/file34ac66f4ddd2
11:57:47 Searching Annoy index using 1 thread, search_k = 3000
11:57:48 Annoy recall = 100%
11:57:49 Commencing smooth kNN distance calibration using 1 thread
11:57:51 Initializing from normalized Laplacian + noise
11:57:51 Commencing optimization for 500 epochs, with 144690 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:57:56 Optimization finished

brain <- FindNeighbors(object = brain, reduction = 'lsi', dims = 2:30)

Computing nearest neighbor graph
Computing SNN

brain <- FindClusters(object = brain, verbose = FALSE, algorithm = 3)
DimPlot(object = brain, label = TRUE) + NoLegend()

gene.coords <- ensembldb::genes(EnsDb.Mmusculus.v79, filter = ~ gene_biotype == "protein_coding")
lncRNA.coords <- ensembldb::genes(EnsDb.Mmusculus.v79, filter = ~ gene_biotype == "lincRNA")
gene.coords <- c(gene.coords,lncRNA.coords)

seqlevelsStyle(gene.coords) <- 'UCSC'
genebody.coords <- keepStandardChromosomes(gene.coords, pruning.mode = 'coarse')

# Flatten the overlapping genes and extend by 2kb upstream of promoters
genebody.coords.flat <- GenomicRanges::reduce(x = genebody.coords)
genebodyandpromoter.coords.flat <- Signac::Extend(genebody.coords.flat,upstream = 2000)

# Retrieve gene names from the original annotation (lost because of flatenning)
genebodyandpromoter.coords.flat$name<- gene.coords[nearest(genebodyandpromoter.coords.flat,genebody.coords)]$gene_name

gene.activities <- FeatureMatrix(
  fragments = fragment.path,
  features = genebodyandpromoter.coords.flat,
  cells = colnames(brain),
  chunk = 20
)

Extracting reads overlapping genomic regions

  |                                                  | 0 % ~calculating  
  |+++                                               | 5 % ~17m 06s      

markers       <- read.csv2(file=paste0(params$out_prefix,'Sten_RNA/clustering/sten_RNA_markers.csv'),row.names = 1)
assay='GA'


brain.markers.ls <- lapply(unique(markers$cluster),function(x){
  marker.genes <- markers[markers$cluster==x & markers$p_val_adj < 0.05 & markers$avg_logFC > 0,"gene"]
  marker.genes <- head(marker.genes,200)
  marker.genes.agg.exp <- Matrix::colSums(brain[[assay]]@data[rownames(brain[[assay]]) %in% marker.genes,])
  marker.genes.agg.exp <- marker.genes.agg.exp / Matrix::colSums(brain[[assay]]@counts)
  return(marker.genes.agg.exp)
})

names(brain.markers.ls)  <- gsub("-",".",unique(markers$cluster))
markers.agg.df           <- do.call(cbind,brain.markers.ls)
colnames(markers.agg.df) <- paste0("marker_",gsub(" ",".",colnames(markers.agg.df)))


brain <- AddMetaData(brain,metadata = markers.agg.df,col.name = colnames(markers.agg.df))



lapply(as.character(colnames(markers.agg.df)),function(x){
  FeaturePlot(brain,x,max.cutoff = quantile(markers.agg.df[,x],0.9,na.rm=TRUE),min.cutoff = quantile(markers.agg.df[,x],0.05,na.rm=TRUE)) + scale_color_viridis_c() + NoLegend()
  #FeaturePlot(brain,x,max.cutoff = 0.08,min.cutoff = 0.02) + scale_color_viridis_c(limits=c(0.02,0.08)) 
  })

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

[[20]]

[[21]]

[[22]]

[[23]]

[[24]]

[[25]]

[[26]]

[[27]]

[[28]]

[[29]]

[[30]]

[[31]]

[[32]]

[[33]]

[[34]]

[[35]]

[[36]]

[[37]]

[[38]]

[[39]]

brain <- RenameIdents(brain,'4' = "mOL",'8' = 'mOL','12' = 'OPC','5' = "Astrocytes",'10'='Astrocytes')
DimPlot(brain,label=TRUE) + NoLegend()

rr paste0(params\(out_prefix,params\)antibody,/bigwig/cluster_,cluster,.bw)

[1] \~/mount/CT/git_test/results/scATAC_P50/bigwig/cluster_olig.bw\

f

Error: object 'f' not found

rr saveRDS(object = brain,file = paste0(params\(out_prefix,params\)antibody,/clustering/clustering_scATAC.Rds))